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Abstract 

Let Q be an open, simply connected, and bounded region in R'', d > 2, 
and assume its boundary dQ is smooth. Consider solving the eigenvalue 
problem Lu = Xu for an elliptic partial diflterential operator L over Q 
with zero values for either Dirichlet or Neumann boundary conditions. We 
propose, analyze, and illustrate a 'spectral method' for solving numerically 
such an eigenvalue problem. This is an extension of the methods presented 
earlier in [5], [6]. 

1 INTRODUCTION 

We consider the numerical solution of the eigenvalue problem 
Lu{s) = - V A (a^ {s)^^\ + j{s)u{s) = Xu{s), s^ncW^ (1) 

with the Dirichlet boundary condition 

u{s) EE 0, s e dn. (2) 

Assume d > 2. Let 17 be an open, simply-connected, and bounded region in 
R"*, and assume that its boundary d^l is smooth and sufficiently differentiable. 
Similarly, assume the functions 7(5) and aij{s), 1 < i, j < d, are several times 
continuously differentiable over fl. As usual, assume the matrix A{s) — [aij{s)] 
is symmetric and satisfies the strong ellipticity condition, 

^'^A{s)^>coe^, sen, ^em." (3) 
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with Co > 0. For convenience and without loss of generahty, we assume 7(5) > 0, 
s G fl; for otherwise, we can add a multiple of u{s) to both sides of ([T]), shifting 
the eigenvalues by a known constant. 

In the earlier papers [S] and [5] we introduced a spectral method for the 
numerical solution of elliptic problems over with Dirichlet and Neumann 
boundary conditions, respectively. In the present work, this spectral method 
is extended to the numerical solution of the eigenvalue problem for ((T])-([2]), and 
in a later section it is also extended to the Neumann problem 

— Au(s) = Au(s), s e f2 

1^ = 0, sedn. 
on 

2 The Dirichlet problem 

Our spectral method is based on polynomial approximation on the unit ball 
in M.'^. To transform a problem defined on to an equivalent problem defined 
on Bd, we review some ideas from |5] and [6,, modifying them as appropriate 
for this paper. 

Assume the existence of a function 

$ : izi n (4) 

onto 

with $ a twice-difFerentiable mapping, and let ^ = : i-^ Bd- For 

onto 

vix) = V (x)) , xeBrfCK'* (5) 

and conversely, 

vis) ^ V {s)) , seHcR'^. (6) 
Assuming v e (fi), we can show 

\/^v{x) = J{xf\7sv{s), s = $(x) 

with J (x) the Jacobian matrix for $ over the unit ball Bd, 



J{x) EE (x) 



dxj 



X G Bd. (7) 



To use our method for problems over a region 17, it is necessary to know explicitly 
the functions i> and J. We assume 

det J{x) ^0, a; e 'Bd- (8) 

Similarly, 

Vsv{s) = K{s)'^V^v{x), X = 4'(s) 
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with K{s) the Jacobian matrix for ^' over fl. By differentiating the identity 

* ($ (a;)) = X, xe Bd 

we obtain 

K {x)) = J {x)~'^ . (9) 

Assumptions about the differentiabihty of v [x) can be related back to assump- 
tions on the differentiabihty of v{s) and <i>(x). 

Lemma 1 // $ e C*^ (Ed) and v e C™ (H) , f/ien w € C"? (Bd) with q = 
min{fc, m}. 

Proof. A proof is straightforward using ■ 

A converse statement can be made as regards v, v, and ^' in ([S]). 
Consider now the nonhomogeneous problem Lu = f, 



k,i=l ^ ^ 



+ j{s)u{s) = f{s), sencR'^. (10) 



Using the transformation it is shown in Thm 2] that (fTU]) is equivalent 
to ^ 

i^kA^) det iJ{x)) ) + [7(2^) det J(a;)] u{x) 

^%dxk V dx, ) (11) 

= / (x) det J(a;), x £ Bd 

with the matrix A (x) = [aij{x)] given by 

A{x) = Jix)'^ A{<i>{x))Jix)-^ . (12) 

The matrix A satisfies the analogue of but over Bd- Thus the original 
eigenvalue problem ([l|)-(l2]) can be replaced by 

- E (^kA^) det {Jix)) + [^(a,) det J(a;)] 

= \u{x) det J(a;), x £ Bd 



k,e=i 



As a consequence of this transformation, we can work with an elliptic problem 
defined over Bd rather than over the original region Q. 

2.1 The variational framework 

To develop our numerical method, we need a variational framework for ([1 
with the Dirichlet condition u = on dH,. As usual, multiply both sides of ([1 



3 



by an arbitary v G Hq (fl), integrate over f2, and apply integration by parts. 
This yields the problem of finding u € Hq (fl) such that 

A {u, v) = (/, v) = i {v) , for all v £ (17) (14) 

with 

ds, v,w e Hq (Q) . 

(15) 

The right side of (|14p uses the inner product (•, •) of (fi). The operators L 
and A are related by 

{Lu,v) ^ A{u,v) , ueH^iQ), veH^in), (16) 

an identity we use later. The function A is an inner product and it satisfies 

\A{v,w)\<ca\\v\\,\\w\\,, v.weH'oi^) (17) 

Aiv,v)>ce\\v\\l veH^in) (18) 

for some positive constants ca and Cg. 
Associated with the Dirichlet problem 

Lu{x)^f{x), xen, feL^in) (19) 
u{x) = 0, xedn (20) 

is the Green's function integral operator 

u{x)^gfix). (21) 

Lemma 2 The operator Q is a bounded and self-adjoint operator from (fl) 
into Hq (fl). Moreover, it is a compact operator from (fl) into Hq (il), and 
more particularly, it is a compact operator from Hq (fl) into Hq (il). 

Proof. A proof can be based on [Tni §6.3, Thm. 5] together with the fact that 
the embedding of Hq {fl) into Hq (fl) is compact. The symmetry follows from 
the self-adjointness of the original problem p^ - ipO]) . ■ 

We convert (|16p to 

{f,v)^A{gf,v), veH^{n), feL'in). (22) 

The problem P^ - (|^D|) has the following variational reformulation: find u G 
H^ (n) such that 

A{u,v)=e{v), yveH^iQ). (23) 

This problem can be shown to have a unique solution u by using the Lax- 
Milgram Theorem to imply its existence; see [Tj Thm. 8.3.4]. In addition, 

Mi<^\\e\\ 

with 1 1 ^11 denoting the operator norm for i regarded as a hnear functional on 



A {v, w) 



dv{s) dw{s) 
dse dsk 



+ j{s)v{s)w{s) 
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2.2 The approximation scheme 

Denote by n„ the space of polynomials in d variables that are of degree < n: 
p £ n„ if it has the form 

\i\<n 

with i a multi-integer, i = (ii, . . . , z^), and \i\ = ii + ■ ■ ■ + id- Over Bd, our 
approximation subspace is 

Xn = {(l-\\x\\l)p{x)\pelln} (24) 

with ||a;||2 = + • • • + x^. The subspaces n„ and Xn have dimension 
However our problem (|14p is defined over fi, and thus we use a modification of 

Xn'. 

A'„ = {^(s)=?A(*(s)):^eAi,} (25) 

The finite dimensional set Xn C 77q (51). This set of functions is used in the 
initial definition of our numerical scheme and for its convergence analysis; but 
the simpler space Xn is used in the actual implementation of the method. 

To solve (|23|) (and thus p9|) - ((20| ) approximately, we use the Galerkin method 
with trial space Xn to find u„ G Xn for which 

AiUn,v)^eiv), WveXn. (26) 

For the eigenvalue problem find u„ G Xn for which 

A{un,v) ^ X{un,v) , VweA'^. (27) 

Write 

N 

with {V'jl^i a- basis of Then (P7)) becomes 

W AT 

^aj.4(V'„^A,) = A^a,(V'„^.), i = l,...,7V (29) 

The coefficients can be related back to a polynomial basis for Xn and to 
integrals over Bd- Let denote the basis of Xn corresponding to the basis 
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{ipj} for Xn. Using the transformation s = <I>(x), 



{ipj,^pt) = / iJj (s)'0j (s) ds 
Jn 



ijjj {x) ipi (x) |det J {x)\ dx 



(30) 



k,e=i 



ds 



ds 



[{V,^, {s)f A{s){\/,^, {s)} + j{s)ij,{s)^,{s) 
+7(a;)-0j(x)?/;i(a;) |detJ(x)|c?x 

Va;'0i (a;)"'" A{x)Vx'ifj (x) + l{x)'ipi [x) {x) |det J (a::)| dx 



with the matrix A{x) given in (fT2|) . With these evaluations of the coefficients, it 
is straightforward to show that ([29l) is equivalent to a Galerkin method for (I12p 
using the standard inner product of (Bd) and the approximating subspace 

2.3 Convergence analysis 

The scheme (j29p is implicitly a numerical approximation of the integral equation 
eigenvalue problem 

XQu — u. (31) 

Lemma 3 The numerical method ^27\ l is equivalent to the Galerkin method 
approximation of the integral equation 113 with the Galerkin method based on 
the inner product A {■, •) for {Q). 



Proof. For the Galerkin solution of f3T|) we seek a function u„ in the form 
, and we force the residual to be orthogonal to . This leads to 



N 



N 



X^ajA{gtpj,ipi) = ^ajA{tpj,'ipi) 



(32) 



for I = 1, ... , N. From fF^ . we have A [Qil^jTi^i) = (V'j, ^'^^ ihns. 

N N 



6 



This is exactly the same as (P^ . ■ 

Let Vn be the orthogonal projection of Hi {B) onto A'„, based on the inner 
product A {■,■)■ Then ((32|l is the Galerkin approximation, 



(33) 



for the integral equation eigenvalue problem (jSip . Much is known about such 
schemes, as we discuss below. The conversion of the eigenvalue problem (|27p 
into the equivalent eigenvalue problem (|33p is motivated by a similar idea used 
in Osborn [25j. 

The numerical solution of eigenvalue problems for compact integral operators 
has been studied by many people for over a century. With Galerkin methods, 
we note particularly the early work of Krasnoselskii [20} p. 178]. The book of 
Chatelin [14j presents and summarizes much of the literature on the numerical 
solution of such eigenvalue problems for compact operators. For our work we 
use the results given in [5] , for pointwise convergent operator approximations 
that are collectively compact. 

We begin with some preliminary lemmas. 

Lemma 4 For suitable positive constants ci and C2, 

ci\\v\\hI(b^) < \MHl,{n) < c2\\v\\hI,{b^) 

for all functions v G Hi (fl), with v the corresponding function of Thus, 
for a sequence in Hi (f2), 



Vn ^ V in Hi (fl) <^=^ Vn ^ V in Hi (Bd) 
with {vn} the corresponding sequence in Hi (Bd). 



(34) 



Proof. Begin by noting that there is a 1-1 correspondence between Hi (fi) and 
Hi {Bd) based on using (g])-®. Next, 



Mh. 



\Vvis)\' + \vis)\' 



ds 



Vv {xf J {x)~'^ J {xy^ Vv (x) +\v{x)\^] |detJ(x)| dx 



< 



max |det J{x)\ 



max ■i max II J (x) "'^ |P, 1 
' xeB 



\Vv {x)\' + \v{xy 



(35) 



dx 



l!w||ffi(n) < C2||w||Hi(ij,) 

for a suitable constant C2 (^l). The reverse inequality, with the roles of ||'y||_ff,5(Bd 
and ||t'||iyi(o) reversed, follows by an analogous argument. I 



Lemma 5 The set Un>iX„ is dense in Hi (ft). 
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Proof. The set U„>iA'„ is dense in {Bd), a result shown in [SJ see (15)]. We 
can then use the correspondence between Hq {fl) and Hq (Bd), given in Lemma 
m to show that U„>iA'„ is dense in (fi) . ■ 



Lemma 6 The standard norm \\ ■ ||i on Hq (fi) and the norm \\v\\a = \/ A {v, v) 
are equivalent in the topology they generate. More precisely, 

V^e\\v\\i<\\v\\A<V^\\v\\i, veH^in). (36) 

with the constants ca, Ce taken from |J7| ) and {_?<$|). respectively. Convergence 
of sequences {wn} is equivalent in the two norms. 

Proof. It is immediate from and p7|) . ■ 



Lemma 7 For the orthogonal projection operator Vn, 

'PnV V as n —> oo, for all v e Hi {Q) . (37) 

Proof. This follows from the definition of an orthogonal projection operator 
and using the result that U„>iA'„ is dense in Hq (fl). M 

Corollary 8 For the integral operator Q , 

||(/-'Pn)^|| ^ as n^oo (38) 

using the norm for operators from Hq (il) into Hq (il). 

Proof. Consider Q and Vn as operators on H^ (fl) into Hq (fl). The result 
follows from the compactness of Q and the pointwise convergence in ([37]); see 
(H Lemma 3.1.2]. ■ 



Lemma 9 {VnG} is collectively compact on Hq (fl) . 

Proof. This follows for all such families {VnG} with Q compact on a Banach 
space y and {Vn} pointwise convergent on y. To prove this requires showing 

{VnGv I ||v||i < 1, n> 1} 

has compact closure in Hq ($7). This can be done by showing that the set is 
totally bounded. We omit the details of the proof. ■ 

Summarizing, {VnG} is a collectively compact family that is pointwise con- 
vergent on Hq (fi). With this, the results in [5], [3] can be applied to ([55]) as 
a numerical approximation to the eigenvalue problem (j3ip . We summarize the 
application of those results to ([33]) . 
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Theorem 10 Let A be an eigenvalue for the problem ([2P"(0); say of multiplicity 
V, and let x^^\ • ■ • i x'"^"* ^6 O' basis for the associated eigenfunction subspace. Let 
£ > Q be chosen such that there are no other eigenvalues of ([2P~0i within a 
distance e of X. Let an denote the eigenvalue solutions of \2T^ that are within e 
of X. Then for all sufficiently large n, say n > hq, the sum of the multiplicities 
of the approximating eigenvalues within (T„ equals v . Moreover, 

max |A-A„|<c max || (/ - P„) X*''' l|i (39) 

Let u be an eigenfunction of associated with A. Let yV„ be the direct 

sum of the eigenfunction subspaces associated with the eigenvalues Xn (z c^n, and 
let ^Un \ . . . , itn'^"'! be a basis for yV„. Then there is a sequence 



fe=i 



for which 



lh-w„|li<c max ||(/-p„)x'^-)||i (40) 

l<k<iy 



for some constant c > dependent on X. 

Proof. This is a direct consequence of results in together with the 

compactness of Q on Hq (fl). It also uses the equivalence of norms given in 

(ES). ■ 



The norms || — "Pn) X*"*^^ II i can be bounded using results from Ragozin 
[25] , just as was done in [5]. We begin with the following result from [SHj. The 
corresponding result that is needed with the Neumann problem can be obtained 
from [9]. 

Lemma 11 Assume w e C''^'^ i^d) for some k > 0, and assume w\gg — 0. 
Then there is a polynomial Qn G '^n for which 

\\w ~qn\\^<D {k, d) n-^ (n-i \\w\\^^^^, + oj (w^^+^ll/n)) (41) 

In this, 

|j|<A;+2 

uj{g,5)= sup \g{x)-g{y)\ 

\x-y\<S 



\i\=k+2 
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Theorem 12 Recall the notation and assumptions of Theorem ] 1 (A Assume the 
eigenfunction basis functions x*-*^^ G (jm+2 ^^^-j assume $ € (7™+2 (^Bd), for 
some m > 1. Then 

max |A- A„| O (rj-™) 

Proof. Begin with ([211) -(gOl). To obtain the bounds for || (/ - Vn) u^'''>\\i given 
above using Lemma [TTl refer to the argument given in [S] . ■ 



3 Implementation 

Consider the implementation of the Galerkin method of (P7|) for the eigenvalue 
problem ([Ij. We are to find the function u„ G Xn satisfying (|29|) . To do so, we 
begin by selecting a basis for n„ that is orthonormal in {Bd), denoting it by 
{^1, . . . , ^pn}, with N = Nn = dimn„. Choosing such an orthonormal basis is 
an attempt to have the matrix associated with the left side of the linear system 
in be better conditioned. Next, let 

Ux)= {l-\\x\\l)^,{x), z = l,...,7V„ (42) 
to form a basis for Xn- As in (pS)) . let {ipi, . . . , i^n} be the corresponding basis 

of Xn- 

We seek 

N 

i=i 

Then following the change of variable s — ^ {x), (|29p becomes 

N ^ 



y~^aj / V'i/'j (a;)"^ ^(a;)V'0i (a;) + "y{x)ipj (x) V'i (a;) |det J (x)| dec 



Xy^aj / -(Aj (a;) V'i (a;) |det J (x)| dec, i = 1, . . . , N 
3=1 



(44) 



We need to calculate the orthonormal polynomials and their first partial deriva- 
tives; and we also need to approximate the integrals in the linear system. For 
an introduction to the topic of multivariate orthogonal polynomials, see Dunkl 
and Xu ^15] and Xu ^30j . For multivariate quadrature over the unit ball in W^, 
see Stroud [28l. 
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3.1 The planar case 

The dimension of n„ is 

Nn = i (71 + 1) (n + 2) (45) 

For notation, we replace x with (x, y). How do we choose the orthonormal basis 
{(pe{x, for n„? Unhke the situation for the single variable case, there are 

many possible orthonormal bases over B = D, the unit disk in R"^. We have 
chosen one that is particularly convenient for our computations. These are the 
"ridge polynomials" introduced by Logan and Shepp [22; for solving an image 
reconstruction problem. We summarize here the results needed for our work. 
Let 

v„ = {P e n„ : (P, Q) - vg e n„_i} 

the polynomials of degree n that are orthogonal to all elements of n„_i. Then 
the dimension of Vn is n + 1; moreover, 

n„ = Vo ® Vi ® • • • ® V„ (46) 

It is standard to construct orthonormal bases of each V„ and to then combine 
them to form an orthonormal basis of n„ using the latter decomposition. As 
an orthonormal basis of V„ we use 

_ 1 TT 

(Pn,k{x,y) = —^Un{x cos {kh) + y sin (kh)) , {x,y) e D, h= — — (47) 
y/TT n+1 

for fc = 0, 1, . . . , n. The function C/„ is the Chebyshev polynomial of the second 
kind of degree n: 

Unit) ^ ^^^'^i^+y^ ^ t = cos6i, -1<<<1, 71 = 0,1,... (48) 
smO 

The family {^n.k}2^Q is an orthonormal basis of V„. As a basis of n„, we order 
{^n,k} lexicographically based on the ordering in (|^7)l and 

{'Pejf^l = {^0,0, ^1,0, 'Pi,!, ^2,0, ■ ■ • , (Pn,0, ■ ■ ■ , ^n,n} 



Returning to p2)) . we define 

i^n,k{x,y) = (1 - a;^ - y^) (pn,k{x,y) (49) 

To calculate the first order partial derivatives oi il)n,k{x^y), we need C/„(t). The 
values of J7„(i) and J7„(i) are evaluated using the standard triple recursion 
relations 

Un+l{t)^2tUn{t)-Un-l{t) 

Ui+,{t) = 2C/„(t) + 2tui{t) - Ui_,{t) 
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For the numerical approximation of the integrals in (I44p , which arc over B 
being the unit disk, we use the formula 



j g{x, y)dxdy^Y.Y.g U ^) '^'^n (50) 

Here the numbers uji are the weights of the {q + l)-point Gauss-Legendre quadra- 
ture formula on [0, 1]. Note that 

.1 9 

/ p{x)dx ^^p{ri)uju 

for all single- variable polynomials p{x) with deg [p) <2q + \. The formula (|50p 
uses the trapezoidal rule with 2g-f-l subdivisions for the integration over B in the 
azimuthal variable. This quadrature (|50p is exact for all polynomials g S ^2q- 
This formula is also the basis of the hyperinterpolation formula discussed in 



3.2 The three— dimensional case 

In 'M? , the dimension of n„ is 

Ar„= (^" + 3^ =i(n+l)(n + 2)(n + 3) 
Here we choose orthonormal polynomials on the unit ball as described in [T5], 

^m,j.A^) = CrnjPj (2||.t|| - I) S fj ,m-2j I^J^^ j 

= Cr.,\\x\r-'^pf^"'-'^^'^\2\\xr l)5,,™-2, (^) , (51) 

j = 0,...,Lm/2j, /3 = 0,1,...,2(to-2j), m = 0,l,...,n 

Here Cm.j = 23+^^-' is a constant, and p^*''™ ^j+a)^ j ^ j^^^ ^^.^ ^j^^ normalized 
Jabobi polynomials which arc orthonormal on [— 1 , 1] with respect to the inner 
product 

(w, w)^ J (1 + ty''-^'+h{t)w{t) dt, 

see for example [T], |17| . The functions Sp,m-2j are spherical harmonic func- 
tions, and they are given in spherical coordinates by 



cos(f</>)T^^(cos6'), /3even 
sin(^0)rfc^ (cos 6*), /? odd 
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The constant C/s^k is chosen in such a way that the functions are orthonormal 
on the unit sphere 5*^ in M^: 

^ Sp^kix) S^j^ix) dS = Sf^^^S^j: 

The functions are the associated Legendre polynomials, see [TO], [13] • Ac- 
cording to we define the basis for our space of trial functions by 

V'mj^/jla;) = (1 - \\x\\'^)ip,n,j.J3{x) 

and we can order the basis lexicographically. To calculate all of the above func- 
tions we can use recursive algorithms similar to the one used for the Chebyshev 
polynomials. These algorithms also allow the calculation of the derivatives of 
each of these functions, see [T7], [32] 

For the numerical approximation of the integrals in (|44p we use a quadrature 
formula for the unit ball B 

g{x)dx^ / / / 9{r,0,<p)r^ sm{(f)) d(f>d9 dr Qq[g] 
B Jo Jo Jo 

QqWl '-^ l^l^l^-^j ''kg ( > arccos(O) I (52) 

Here ^(r, 0, 0) = g(x) is the representation of g in spherical coordinates. For the 
d integration we use the trapezoidal rule, because the function is 27r— periodic 
in Q. For the r direction we use the transformation 

^'"^'^ = L (^) " (^) ? 

' 1 , fCk + i 



k=l^ 



where the v'j. and (^k are the weights and the nodes of the Gauss quadrature 
with q nodes on [— 1 , 1] with respect to the inner product 

{v,w)=J {l + tfv{t)w{t)dt 

The weights and nodes also depend on q but we omit this index. For the </> 
direction we use the transformation 

/•I 

J-i 
1 



ujjv(arccos{£,j)) 
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Figure 1: Images of ((55)) . with a = 0.5, for lines of constant radius and constant 
azimuth on the unit disk. 



where the tuj and are the nodes and weights for the Gauss-Legendre quadra- 
ture on [—1,1]. For more information on this quadrature rule on the unit ball 
in K^, see [28] , 

Finally we need the gradient in Cartesian coordinates to approximate the 
integral in (|44p . but the function ipm,j.p{x) in (|5ip is given in spherical coordi- 
nates. Here we simply use the chain rule, with x — [x, y, z), 

d d d sin(^?) 

■^v{r, 0, 4>) = T-v{r, 6, 0) cos(6') sm{(f>) - T^w(r, P 



dx ' ' 9r ' ' 09 ' ' rsin((/)) 

^ f a ,nCOs(6I)cos(0) 



and similarly for ^ and ^ 



4 Numerical example 

Our programs are written in Matlab. The transformations have been so chosen 
that we can invert explicitly the mapping <I>, to be able to better construct our 
test examples. This is not needed when applying the method; but it simplified 
the construction of our test cases. The eigenvalue problem being solved is 

Lu(s) = -Au = Aii(s), sencR'^ (53) 
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t 



-2 -2 



s 



Figure 2: Eigcnfunction corresponding to the approximate eigenvalue 
Ai = 2.96185. 



which corresponds to choosing A = I. Then we need to calculate 

I(x) = J(x)-i J(x)-^ 

4.1 The planar case 



(54) 



For our variables, we replace a point x Cz B2 with {x, y), and we replace a point 
s e r2 with (s, t). Define the mapping ^ : B2 hy {s,t) — ^ {x, y), 

(55) 



s = X — y + ax 
t — X + y 



with < a < 1. It can be shown that $ is a 1-1 mapping from the unit disk B. 
In particular, the inverse mapping ^' : il — )■ i? is given by 

1 r 



a 

1 r 



y 



1 + + a (5 + t) 
at - (^-1 + ^/l + a(s + i)) 



(56) 



In Figure [TJ we give the images in il of the circles r — j /lO, j — 1, . . . ,10 and 
the azimuthal lines 9 = j'tt/IO, j — 1, . . . , 20. 

The following information is needed when implementing the transformation 
from — Au — Xu on f2 to a new equation on B2: 



l + 2ax -1 
1 1 
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Figure 3: Eigcnfunction corresponding to the approximate eigenvalue 
A2 = 7.24761. 



det (J) -2(1 



A = j{x) ^j{xy 



2(1 



ax) 
1 



1 

-1 



1 

1 + 2ax 



1 

ax 



2a^x^ + 2ax + 1 



2(1 + ax)' 

We give an example for this region O with a — 0.5. Figures [2] and [3] contain 
the computed eigenfunctions for the two smallest eigenvalues; these are based 
on the degree n = 8 approximation. 

Because the true eigenfunctions and eigenvalues are unknown for almost all 
cases (with the unit ball as an exception), we used other methods for studying 

experimentally the rate of convergence. Let aI'^^ denote the value of the fc*'' 
eigenvalue based on the degree n polynomial approximation, with the eigenval- 
ues taken in increasing order. Let vln^ denote a corresponding eigcnfunction, 

SfW (x)=^a5"V,(x) 



with = 



the eigenvector of (|44)) associated with the eigen- 



value A„ . We normalize the eigenvectors by requiring ||a^"^||oo = 1- Define 



A. 
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O Eigenvalue A,'^' 
Eigenvalue A,'^' 




Figure 4: The values ol 



,(fc) ,(fc) 



for A; = 1, 2 for increasing degree n. 



Figures [4] and [5] show the decrease, respectively, of A„ and £)„ as n increases. 
In both cases, we use a semi-log scale. Also, consider the residual 

Figure [5] shows the decrease of ||i?l'^'*||oo, again on a semi-log scale. 

These numerical results all indicate an exponential rate of convergence as a 
function of the degree n of the approximations |Arf' : ^ > l| and |urf' : n> l|. 

In Figure m the maximum accuracy for A^^-* appears to have been found with 
the degree n = 12, approximately. For larger degrees, rounding errors dominate. 
We also see that the accuracy for the first eigenvalue-eigenfunction pair is better 
than that for the second such pair. 

4.2 The three— dimensional case 

Here we consider the problem of finding eigenvalues and eigenfunctions for the 
Neumann problem in 17 C M-^: 

— Am(s) = Am(s), s G il 
on 
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O Eigenfunction u*^' 
^ Eigenfunction u*^' 




Figure 5: The values of H'ulj^ij'.i — ui^^ ||oo for fc = 1, 2 for increasing degree n. 



Problem ([57]) is equivalent to 

-Au{s)+u{s) = (A + 1)m(s), sei} 

0, sedn 



du{s) 



(58) 



dn 

and + / : i-^ L^(ri) is an invertible self-adjoint operator with 



Dr 



du{s) 
dn 



0,sedn 



So there is a continuous solution operator G : L^(fi) i-^ Dl, such that 

(-A + /)oG|l2(j,) =/ 

with I the identity operator on (fl). If we consider G : H^{n) i-^ H^{U), then 
G is a compact operator, because of the compact imbedding H^{n) ^ L^{^1) 
or ^ see |2T] or 

We follow now Section 12.11 to present the variational framework. A solution 
of the inhomogeneous problem 



(59) 



satisfies 



I \ ^uis) + uis)] vis) ds ^ [ f(s)v(s)ds toi all v e H\n) 
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Figure 6: The values of ||i?n ||oo for k ~ 1,2 for increasing degree n. 

Applying integration by parts and using the fact that the normal derivative of 
u ^ Dl is zero on d^l we derive 

/ \I su{s)\7 sv{s) + u{s)v{s) ds^ I f{s)v{s) ds for all v £ H\n) 
Jn Ja 

We denote the left hand side of this equation by A{u, v) and from the Cauchy- 
Schwartz inequality we derive 

A{u,v) < \\u\\Hnn) \Hm{n) 

and we have the equality 

Aiu,u) = 

Because we assumed that the boundary dft is at least C^, regularity theory 
shows that a solution u G H^{il) of the variational problem 

A{u, v) = if, v) for aU v e (n) (60) 

fulfills u £ Dl; see again [3T] or [5^. So the problems and (|SD|) are 
equivalent. 

Instead of (|58p we consider the equivalent variational problem to find u G 
iJ^(ri) which solves 

A{u,v)^{X + l) [ u{s)v{s)ds for all u e 
Ja 
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and this is equivalent to 

/ Vsu{s)Vsv{s) ds^X I u{s)v{s) ds for all v £ H^{n) (61) 
Jn Jn 

Equation (j6ip is the starting point for our numerical approximation scheme, see 
also 1271 First we transfer equation (|61[) to an equation on the domain B3 with 
the help of a transformation $ : ^3 51. So (pTj) becomes 



V xu{x)A{x)V xvix)\ det(J(a;)| dx 

(62) 

= A / u{x)v{x)\ det(J(a;)| dx for all v € ^^(^3) 

where A{x) — J {x)~^ J {x)~^ ] see ©-([T]) for the definition of the functions and 
J{x). According to Section we need a sequence of subspaces C X^+i C 
iJi(B3) with 

00 

n=l 

Because there are no boundary conditions imposed on H^{B^) we can use 

Xn = {p{x) I p e n„} 

where n„ is the space of polynomials in 3 variables of degree n or less. As a 
basis we choose 



{ip,{x)\i = l,...,Nn}, Nn 



n + 3 
3 



where {(pi} is an enumeration of the orthogonal basis {(pm,j,i3} given in (|5ip. To 
approximate the solutions u{x) of we use 

«l*Ha;)=^af ^,(x),:z =,..., ^„ 

and the coefficients aj*'' for the eigenvalue approximation A^*-* are given as solu- 
tions of the finite eigenvalue problem 



( / Vx^jix)A{x)Vxfkix)\det{J{x))\dx 



a ■ 



¥3j(x)¥3fc(x)| det(J(a;))| dx I a'-''\ k ^ 1, . . . , N„ 
S3 / 

The functions ^^^jix) can be calculated explicitly and all integrals in formula 
([63)) are approximated by the quadrature formula ((52)) with q = n. The conver- 
gence analysis of Section [2.31 can be used without any modifications. 
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Figure 7: The boundary of fli 



To test our method we use two different domains. Let B3 denote the closed 
unit ball in R'^. The domain fli = ^1(33) is given by 



so i?3 is transformed to an ellipsoid ili; see Figure [T] The domain fl2 is given 



where we used polar coordinates {p, <j), 9) E [0, 1] x [0, 2tt] x [0, tt] to define the 
mapping $2. Here the function S : S'^ = dB-s 1-^ (l,oo) is a function which 
determines the boundary of a star shaped domain The restriction S{4i, 6) > 
1 guarantees that $2 is injective, and this can always be assumed after a suitable 
scaling of CI2 ■ For our numerical example we use 




Xi — 3X2 
2xi + X2 
X1+X2+ x^ 



by 




(64) 




Finally the function t is defined by 
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Tabic 1: Numerical results for fli, h = 0.0001 to approximate R, 



71 










■^V 7 15 / 


^'71 




1 


4 


4.26^; - 2 


i.oo£;- 1 


9.93E - 2 


1.31^-1 


1A5E- 1 


2.63^- 1 


2 


10 


4.26i; - 2 


i.ooe;- 1 


9.93£; - 2 


1.31^;- 1 


1A5E-1 


2.63^; - 1 


3 


20 


1.42£;-4 


5.67S-4 


5.47S - 3 


l.OlS-2 


2.28S- 2 


5.15^-2 


4 


35 


1.42^; - 4 


5.67£; - 4 


5.47£; - 3 


1.01£;-2 


2.28£;- 2 


5.15^;- 2 


5 


56 


1.06i;- 7 


8.38£;- 7 


1.04£;-4 


2.72^; - 4 


1.22£;- 3 


3.54E - 3 


6 


84 


1.06^-7 


8.38S- 7 


1.04S-4 


2.72S-4 


1.22^-3 


3.54E - 3 


7 


120 


2.53E' - 11 


4.31E' - 10 


1.24^-6 


4.85S-6 


3.02^-5 


1.08^;- 4 


8 


165 


2.53E;- 11 


4.31^;- 10 


1.24^;- 6 


4.85^;- 6 


3.02^;- 5 


1.08£;- 4 


9 


220 


2.22i;- 11 


6.78E - 14 





6.32£; - 8 


4.25£;- 7 


1.80i;-6 


10 


286 


4.47£;- 11 


1.81£;- 13 





5.77^-8 


4.21S- 7 


1.80S- 6 


11 


364 


1.84£; - 13 


5.19£;- 13 








1.48^;- 8 


4.15^-8 


12 


455 


2.07i; - 13 


l.lSi;- 13 








2.21E-Q 


1.88^;- 8 


13 


560 


1.52i;- 13 


i.giE- 13 








5.81S- 9 


3.43S- 8 


14 


680 


4.64^; - 13 


5.56^; - 14 








1.21^-8 


4.26^;- 8 



Table 2: Numerical results for Q.2 



n 


Nn 


An ' — \\^\ 


lA^^-Al^^l 






1 


4 


i.mE - 1 


3.21£;- 1 


2ME- 1 


3.31S- 1 


2 


10 


3.60^; - 1 


3.21£;- 1 


2.86£; - 1 


3.31£;- 1 


3 


20 


8.16i;- 2 


8.53£;- 2 


8.32£; - 2 


1.05^;- 1 


4 


35 


8.16£;- 2 


8.53S- 2 


8.32S- 2 


1.05S- 1 


5 


56 


1.99£; - 2 


2.27£; - 2 


2.84£; - 2 


3.11£;-2 


6 


84 


1.99i;- 2 


2.27£;- 2 


2.84£; - 2 


3.11^;- 2 


7 


120 


1.48£;- 2 


1.56S- 2 


2.49S- 2 


2.71£:- 2 


8 


165 


1.48^-2 


1.56S- 2 


2.49S- 2 


2.71S-2 


9 


220 


4.77E; - 3 


5.96^;- 3 


1.14^;- 2 


1.47^;- 2 


10 


286 


4.77i; - 3 


5.96^; - 3 


1.14^;- 2 


1.47^;- 2 


11 


364 


8.34£; - 4 


1.28S- 3 


3.25B- 3 


4.76S-3 


12 


455 


8.34E' - 4 


1.28^;- 3 


3.25S - 3 


4.76_B-3 


13 


560 


1.88£: - 4 


2.55£: - 4 


1.33£: - 3 


1,G2£: - 3 


14 


680 


1.88£; - 4 


2.55^; - 4 


1.33^; - 3 


1.62£;-3 
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Figure 8: A view of dfl2 



where the exponent 5 impUes $2 G C^(Bi{0)). See [6J for a more detailed 
description of $2; one perspective of the surface is shown in Figure [8] 

For each domain we calculate the approximate eigenvalues A^*^ , Xn^ = < 
An ^ < < ■ ■ ■ and eigenfunctions Un\ i = 1, . . . , Nn, for the degrees n = 
1, . . . , 15 (here we do not indicate dependence on the domain fl). To analyze 
the convergence we calculate several numbers. First we estimate the speed of 
convergence for the first two eigenvalues by calculating |A^g — Xn^\, i = 1,2, 
n— 1, . . . , 14. Then to estimate the speed of convergence of the eigenfunctions 
we calculate the angle (in L^{n)) between the current approximation and the 
most accurate approximation /(un^w^g), i = 1,2, n — 1, . . . , 14. Finally, an 
independent estimate of the quality of our approximation is given by 

i?W^|-AuW(s)-A««(s)| 

where we use only one s G fi, given by <i>(l/10, 1/10, 1/10). To approximate 
the Laplace operator we use a second order difference scheme with h — 0.0001 
for ill and h = 0.01 for $12. The reason for the latter choice of h is that our 
approximations for the eigenfunctions on CI2 are only accurate up three to four 
digits, so if we divide by the discretization errors are magnified to the order 
of 1. 

The numerical results for fli are given in table 14.21 The graphs in Figures 
[QHTTI seem to indicate exponential convergence. For the graphs of Z{un\u^il), 
see Figure [TOl We remark that we use the function arccos(a:) to calculate the 
angle, and for n w 9 the numerical calculations give a; = 1, so the calculated 
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Figure 9: f2i: errors 



for the calculation of the first two eigenvalues 



angle becomes 0. For the approximation of i?„ one has to remember that we 
use a difference method of order 0{h?) to approximate the Laplace operator, so 
we can not expect any result better than 10^* if we use h = 0.0001. 

As we expect, the approximations for ^2 with the transformation <I>2 present 
a bigger problem for our method. Still from the graphs in Figure [T^ and [T51 we 
might infer that the convergence is exponential, but with a smaller exponent 
than for Vli. Because $2 G C^^B^) we know that the transformed eigenfunctions 
on are in general only C^, so we can only expect a convergence of 0(n~^). 
The values of n which we use are too small to show what we believe is the true 
behavior of the Rn\ although the values for n = 10 ... 14 seem to indicate some 
convergence of the type we would expect. 

The poorer convergence for U,2 as compared to VLi illustrates a general prob- 
lem. When defining a surface dVt by giving it as the image of a 1-1 mapping from 
the unit sphere S"^ into R"^, how does one extend it to a smooth mapping from 
the unit ball to fi? The mapping in (j64p is smooth, but it has large changes in 
its derivatives, and this affects the rate of convergence of our spectral method. 
We are working at present on this problem, developing a numerical method to 
find a well-behaved polynomial mapping $ when given only its restriction to 
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-& — Eigenfunction for X' 




Figure 10: f2i: angles Z yun , Ui^j between the approximate eigenfunction u' 
and our most accurate approximation u^^ ~ u^^\ 
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